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FIGURES 


1 .  Plot  showing  the  measured  NML  transmitter  VLF  electric  field  measured  by 
IMAGE  (red)  and  the  end-to-end  model  predicted  fields  (black).  The  quantitative 
agreement  is  excellent  and,  with  our  simulations  of  ionospheric  penetration,  no 
-20-dB  empirical  correction  is  required. 
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2.  Test  of  FDTD  model  convergence  from  high-latitude  NML  transmitter.  The  panels 
show  the  spectrum  of  low-altitude  (top)  and  high-altitude  (bottom)  electric  fields  as 
the  grid  spacing  in  the  simulation  is  reduced  by  25%.  The  change  in  the  low- 
altitude  fields  is  almost  negligible  across  the  entire  band,  indicating  that  the  results 
are  accurate.  The  high-altitude  fields  differ  somewhat  above  20  kHz,  but  the  change 
is  still  small. 
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3.  Test  of  FDTD  model  convergence  from  low-latitude  NPM  transmitter.  The  panels 
show  the  spectrum  of  low-altitude  (top)  and  high-altitude  (bottom)  electric  fields 

as  the  grid  spacing  in  the  simulation  is  reduced  by  25%.  The  change  in  the  low- 
altitude  fields  is  almost  negligible  across  the  entire  band,  indicating  that  the  results 
are  accurate.  However,  the  high-altitude  fields  differ  somewhat  above  5  kHz  and 
differ  enormously  above  15  kHz. 
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4.  The  low-altitude  (90  km)  field  spectrum  for  a  magnetic  latitude  corresponding 
to  the  NML  transmitter  (high  latitude)  for  3  different  grid  spacings. 
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5.  The  high-altitude  (130  km)  field  spectrum  for  a  magnetic  latitude  corresponding 
to  the  NML  transmitter  (high  latitude)  for  3  different  grid  spacings. 
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6.  The  low-altitude  (90  km)  field  spectrum  for  a  magnetic  latitude  corresponding 
to  the  NPM  transmitter  (low  latitude)  for  3  different  grid  spacings. 
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7.  The  high-altitude  (130  km)  field  spectrum  for  a  magnetic  latitude  corresponding 
to  the  NPM  transmitter  (low  latitude)  for  3  different  grid  spacings. 

14 

8.  Low-altitude  fields  produced  by  a  10-kHz  source  computed  using  the  FD  and 
TD  codes.  The  agreement  is  excellent,  validating  the  new  FD  code. 
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9.  High-altitude  fields  produced  by  a  10-kHz  source  computed  using  the  FD  and 

TD  codes.  The  agreement  is  again  excellent.  17 
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10.  Low-altitude  fields  produced  by  a  20-kHz  source  computed  using  the  FD  and 
TD  codes. 
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1 1 .  High-altitude  fields  produced  by  a  20-kHz  source  computed  using  the  FD  and 
TD  codes. 
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1.  SUMMARY 


The  main  objective  of  this  project  is  to  apply  a  full-wave  electromagnetic  field  model 
to  predict  the  injection  of  wave  energy  produced  by  a  low-altitude  very  low  frequency 
(VLF)  transmitter  through  an  arbitrary  and,  therefore,  realistic  ionosphere.  Our  specific 
task  is  to  deliver  an  accurate  model  of  high-altitude  electromagnetic  fields  and  wave 
energy  injection  to  be  used  in  analyzing  and  predicting  controlled  radiation  belt  particle 
loss  to  reduce  potential  threats  to  satellite  systems.  Accurate  prediction  of  high-altitude 
wave  energy  injection  by  low-altitude  VLF  wave  transmission  is  a  challenging  but 
essential  step  towards  accurate  modeling  of  the  impact  of  man-made  VLF  transmissions 
on  natural  radiation  belt  dynamics.  The  technical  challenge  lies  in  simulating  propagation 
in  a  lossy,  anisotropic  and  arbitrarily  inhomogeneous  ionosphere,  for  which  approximate 
techniques  such  as  ray  tracing  do  not  apply. 

Our  technical  approach  is  to  develop  a  full-wave  numerical  model  to  simulate  wave 
power  injected  into  die  ionosphere  from  a  ground-level  VLF  transmitter.  The  main 
advantage  of  a  finite  difference-based  approach  is  its  ability  to  compute  high  altitude 
fields  in  arbitrary  ionospheres,  including  the  influence  of  Earth’s  magnetic  field,  with  no 
implicit  approximations. 

Our  progress  in  the  second  year  of  this  project  was  significantly  hindered  by  delays  in 
receiving  funding  to  support  the  work.  These  funding  delays  also  meant  that  extra  time 
and  effort  was  expended  ramping  the  effort  down  and  back  up  again.  Nevertheless, 
substantial  progress  was  made.  A  complete  end-to-end  run  of  the  high-altitude  VLF 
power  prediction  model  has  been  completed  using  our  modeled  130- km  altitude  VLF 
power  predictions  as  the  input  to  the  higher  altitude  ray-tracing  code.  This  resulted  in 
good  agreement  with  high-altitude  VLF  field  measurements  from  the  IMAGE  satellite. 
Previous  end-to-end  runs  like  this,  which  used  semi-empirical  ionospheric  penetration, 
had  revealed  a  ~20-dB  discrepancy  between  the  measured  and  predicted  fields.  Our 
model  does  not  require  any  such  correction  factor  and  appears  to  show  that  this 
discrepancy  originates  in  file  ionospheric  penetration  calculation. 

We  have  also  made  significant  progress  in  understanding  the  convergence  of  our 
simulations  and,  thus,  also  quantitatively  bounding  the  absolute  errors  in  the  results. 

We  have  run  a  series  of  simulations  using  parameters  corresponding  to  the  NML  (high- 
latitude)  and  NPM  (low-latitude)  transmitters  with  grid  spacings  of  1000,  500,  and  250 
meters.  Finite  difference  simulations  are  guaranteed  to  eventually  converge  to  the  correct 
answer  as  this  grid  spacing  is  reduced,  but  it  is  not  always  easy  to  know  what  grid 
spacing  is  required  to  maintain  errors  below  a  specified  level.  We  have  found  that,  for  a 
given  uniform  grid  spacing,  the  low-altitude  fields  can  be  computed  correctly  while  the 
high-altitude  fields  are  dramatically  incorrect.  This  shows  that  a  nonuniform  grid 
approach,  in  which  low  altitudes  are  resolved  coarsely  and  only  the  highest  altitudes  are 
resolved  finely,  is  the  key  to  achieving  efficient  simulations  of  high-altitude  fields.  At 
high  latitudes  (i.e.,  NML),  500-m  grid  spacing  results  in  fairly  accurate  high-altitude 
fields  in  the  20-25  kHz  band.  But  at  low  latitudes  (i.e.,  NPM),  250  m  and  smaller  grid 
spacing  is  required  to  achieve  the  same  accuracy.  The  convergence  in  all  cases  shows 
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second-order  accuracy,  which  is  expected.  This  also  means  that  correct  answers  can  be 
estimated  by  extrapolating  the  convergence  to  its  limit. 

Lastly,  the  realization  that  very  fine  grid  spacings  are  required  to  achieve  accurate 
results  at  low  latitudes  has  led  us  to  develop  a  more  efficient  frequency  domain  (i.e., 
single  frequency)  VLF  solver.  Significant  progress  has  been  made  in  developing  a  more 
efficient  frequency  domain  solver  and  validating  its  accuracy.  We  have  found,  and  show 
below,  that  our  preliminary  frequency  domain  code  gives  answers  in  quantitative 
agreement  with  the  time  domain  code  but  is  about  90  times  faster  for  single-frequency 
computations.  This  major  enhancement  in  computational  speed  justifies  the  further 
development  needed  for  this  solver. 

The  ongoing  code  development  will  enable  accurate  predictions  of  the  ionospheric 
wave  fields  and  their  statistical  uncertainty.  The  model  output  will  be  designed  for  easy 
integration  into  ray-tracing  codes  that  complete  the  overall  simulation  of  wave  fields  in 
the  radiation  belts.  The  end  result  will  be  the  capability  to  predict  more  accurately  the 
injection  of  VLF  wave  energy  into  the  magnetosphere,  which  is  a  key  component  in 
understanding  the  influence  of  VLF  on  radiation  belt  dynamics. 

2.  BACKGROUND  AND  ACCOMPLISHMENTS 

Predicting  the  wave  power  injected  into  the  magnetosphere  from  a  ground-level  VLF 
or  lower  frequency  transmitter  is  an  essential  component  in  modeling  radiation  belt 
dynamics.  But  accurately  computing  VLF  power  penetration  through  the  ionosphere  is 
not  easy.  Ray  tracing,  an  effective  technique  in  the  magnetosphere,  is  not  a  good 
approximation  because  the  ionosphere  changes  significantly  on  the  scale  of  a  VLF 
wavelength.  Consequently,  wave  reflection,  mode  conversion,  and  other  full-wave  effects 
that  ray  tracing  generally  neglects  are  important.  Full-wave  mode  theory  techniques  [e.g., 
Papperi  and  Ferguson,  1986]  are  effective  for  reliably  computing  subionospheric  fields 
produced  by  VLF  transmitters.  But  mode  theory  techniques  are  difficult  to  apply  above 
approximately  60-70  km  for  realistic  ionospheres  in  which  the  parameters  are  not 
exponentially  varying.  This  upper  altitude  is  too  low  to  be  of  direct  use  to  the  radiation 
belt  problem.  Past  numerical  calculations  have  been  summarized  in  forms  potentially 
useful  for  this  application,  but  these  calculations  have  always  made  assumptions  that 
limit  their  applicability. 

An  effective  and  efficient  approach  to  computing  ionospheric  VLF  fields  under  very 
general  conditions  is  full-wave  finite  difference  simulation.  This  approach  contains 
essentially  no  approximations  other  than  approximating  derivatives  as  discrete 
differences.  Finite  difference  techniques,  in  general,  have  found  wide  application  in 
numerical  simulations  of  scattering,  propagation,  and  other  electromagnetic  phenomena 
[Taflove  and Hagness,  2000],  The  finite  difference  technique  is  well  suited  to  specific 
aspects  of  the  VLF  problem,  namely,  the  relatively  long  wavelengths  involved  and  a 
domain  with  complex  inhomogeneities.  In  fact,  it  is  really  the  only  technique  capable  of 
handling  the  almost  completely  arbitrary  inhomogeneities  that  appear  in  the  VLF 
ionospheric  penetration  problem. 

Its  chief  disadvantage  is  computational  cost;  by  explicitly  computing  the 
electromagnetic  fields  everywhere  in  the  domain  of  interest,  it  is  a  brute-force  approach. 
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Modem  computer  speed,  however,  is  sufficient  that  speed  is  no  longer  a  major  limitation 
for  this  class  of  problem.  And,  most  importantly,  we  have  already  developed  and 
validated  a  code  that  can  compute  VLF  fields  at  ionospheric  altitudes  from  an  essentially 
arbitrary  source  and  in  an  arbitrary  ionosphere  [Cummer,  2000;  Hu  and  Cummer,  2006]. 
We  have  applied  this  code  to  short-  and  long-distance  VLF  propagation  for  a  variety  of 
ionospheric  and  lightning  remote  sensing  applications. 

The  main  objective  of  this  work  is  to  deliver  a  full-wave  electromagnetic  field  model 
to  predict  the  injection  of  wave  energy  produced  by  a  low-altitude  VLF  transmitter 
through  the  ionosphere.  Our  goal  is  to  accurately  model  high-altitude  fields  to  better 
understand  radiation  belt  loss  mechanisms  and  to  determine  whether  they  can  be 
modified  by  VLF  wave  energy  injection  at  key  locations.  The  work  completed  during  the 
second  year  is  described  below  in  Section  4. 


3.  TECHNICAL  APPROACH 


Terrestrial  VLF  transmitters  play  a  significant  role  in  the  depopulation  of  energetic 
particles  from  the  radiation  belts.  Computing  VLF  fields  in  and  above  the  ionosphere 
from  ground  transmitters  is  difficult  because  wave  propagations  in  a  lossy,  anisotropic 
and  arbitrarily  inhomogeneous  ionosphere  imply  that  any  approximate  techniques  such  as 
ray  tracing  and  mode  theory  for  wave  simulation  are  not  accurate.  Our  technique 
approach  is  to  modify  and  apply  an  existing  finite  difference  code  to  this  problem.  This 
code,  by  approximating  derivatives  as  discrete  differences,  solves  everywhere  in  the 
computational  domain  the  Maxwell  equations  coupled  to  the  equation  for  vector  electric 
current  in  a  magnetized  cold  plasma 


di  Q  , 

-£■ +  ^  x  b  b)+s0(OpE, 


where  0)p  is  the  local  plasma  frequency,  ®b  is  the  local  gyrofrequency,  v  is  the  local 

collision  frequency,  and  bs  is  the  unit  vector  pointing  in  the  direction  of  Earth’s  magnetic 
field  [Budden,  1985].  All  of  cold  plasma  magnetoionic  theory  is  contained  in  these 
equations  and,  thus,  all  relevant  physics  of  linear  electromagnetic  waves  are 
automatically  included  in  this  simulation.  For  example,  the  complicated  Appleton- 
Hartree  refractive  index  formula  is  derived  directly  from  the  time  harmonic  equivalent  of 
the  above  equation.  Both  electron  and  ion  effects  are  included  in  our  model  by  separately 
computing  the  current  from  each  particle  species  through  equations  of  the  same  form. 

3.1.  2D  Cylindrical  Symmetric  FDTD  Modeling 

Our  primary  computational  code  remains  a  2D  azimuthally  symmetric  cylindrical 
coordinates  code.  This  FDTD  code  has  been  developed  and  carefully  validated  for  years 
through  comparisons  with  ground  measurements  and  mode  theory  based  LWPC 
simulations  [Hu  and  Cummer,  2006],  There  are  several  challenges  in  developing  a 
functional  code.  The  strong  anisotropic  nature  of  the  medium  made  simply  developing  a 
stable  finite  difference  scheme  difficult.  Also,  the  wavelength  scales  are  strongly  variable 
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across  the  computational  domain  because  the  refractive  index  of  the  escaping  whistler¬ 
mode  wave  energy  is  very  high,  leading  to  very  short  wavelengths  at  altitude  above  100 


km. 


3.2.  Time  Harmonic  Modeling 

For  many  of  the  computations  needed  for  this  project,  only  a  single  frequency  of  a 
small  set  of  frequencies  needs  to  be  simulated.  A  full  time  domain  simulation,  which 
essentially  computes  all  frequencies  simultaneously,  may  be  less  efficient.  Because  of 
the  multiscale  nature  of  the  problem,  in  which  wavelengths  are  much  shorter  in  the 
ionosphere  than  in  the  lower  altitude  free  space  region,  we  would  like  to  have  variable 
grid  spacing  throughout  the  computational  domain.  This  is  difficult  to  implement  in  a 
time  domain  code,  but  it  is  straightforward  to  implement  in  a  frequency  domain  or  time- 
harmonic  code. 

A  time-harmonic  code  is  based  on  a  frequency  domain  representation  of  the  Maxwell 
equations  and  governing  equations  for  a  cold  plasma,  in  which  the  time  derivatives  are 
converted  to  complex  multiplications  by  ja>.  The  equations  are  approximated  by  spatial 
finite  differences,  and  the  result  is  a  time-static,  matrix  equation  for  all  field  components 
at  all  grid  points  in  the  simulation.  Although  the  number  of  unknowns  will  be  too  large 
for  a  direct  matrix  solution,  the  resulting  sparse  matrix  equation  can  be  solved  using  any 
one  of  a  number  of  classes  of  iterative  sparse  matrix  solvers. 

4.  PRIMARY  RESULTS 


Our  technical  accomplishments  from  the  second  year  of  our  research  effort  are 
summarized  below. 

4.1.  End-to-End  Model  Validation 

We  ran  a  careful  simulation  of  the  high-altitude  (130  km)  fields  in  a  2000-km  radius 
circle  around  the  25.2-kHz  NML  transmitter.  Our  AFRL  colleagues  used  these  high- 
altitude  fields  as  input  to  a  ray-tracing  code  that  computes  the  VLF  power  levels  at  higher 
altitudes  in  the  magnetosphere.  The  computed  magnetospheric  electric  field  levels  were 
compared  to  fields  measured  by  the  RPI  instrument  on  the  IMAGE  satellite.  Figure  1 
shows  an  example  of  this  comparison.  The  predicted  and  measured  field  levels  are  in 
good  quantitative  agreement  with  no  significant  bias  towards  higher  or  lower  fields. 
Similarly  good  agreement  was  achieved  with  data  from  several  other  IMAGE  passes. 

These  results  appear  to  confirm  that  the  ~20-dB  correction  factor  needed  when 
approximate  ionospheric  penetration  calculations  are  used  originates  in  the  ionospheric 
penetration  calculation  itself.  Effort  is  ongoing  to  confirm  this  with  measurements  from 
other  satellites,  which  would  eliminate  uncertainty  in  the  measurement  calibration  as  a 
possible  source  of  the  discrepancy.  If  it  holds  up,  this  result  means  that  prior  radiation 
belt  loss  calculations  may  need  to  be  revisited  in  light  of  weaker  ionospheric  penetration 
than  that  predicted  by  simplified  models. 
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Figure  1.  Plot  showing  the  measured  NML  transmitter  VLF  electric  field  measured 
by  IMAGE  (red)  and  the  end-to-end  model  predicted  fields  (black).  The 
quantitative  agreement  is  excellent  and,  with  our  simulations  of  ionospheric 
penetration,  no  ~20-dB  empirical  correction  is  required. 

4.2.  High-Latitude  vs.  Low-Latitude  Model  Convergence — Initial 
Results 

As  noted  in  previous  reports,  simulations  are  indicating  that  there  is  a  significant 
latitude  dependence  on  the  numerical  parameters  needed  for  accurate  simulations.  This  is 
demonstrated  in  Figures  2  and  3,  which  show  initial  results  in  this  effort  for  field 
computations  at  the  NML  and  NPM  transmitter  latitudes,  respectively.  For  both 
latitudes,  we  see  only  minor  differences  in  the  low-altitude  fields  vs.  frequency  as  the 
grid  size  is  reduced.  This  confirms  that  the  low-altitude  fields  are  accurate  across  the 
entire  VLF  band  for  1-km  grid  spacing. 

However,  the  high-altitude  fields  are  different.  The  high-latitude  (NML)  fields 
have  come  close  to  converging  for  0.75-km  grid  spacing,  although  a  finer  grid  may  be 
needed  for  accurate  fields  above  20  kHz.  The  low-latitude  (NPM)  fields,  however,  are 
not  even  close  to  converged.  A  modest  change  in  the  grid  spacing  (1  to  0.75  km) 
increases  the  fields  by  more  than  a  factor  of  10  above  20  kHz.  Clearly,  a  finer  grid 
spacing  is  needed  for  accurate  simulations  at  NPM  latitudes. 
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Figure  2.  Test  of  FDTD  model  convergence  from  high-latitude  NML  transmitter. 
The  panels  show  the  spectrum  of  low-altitude  (top)  and  high-altitude  (bottom) 
electric  fields  as  the  grid  spacing  in  the  simulation  is  reduced  by  25%.  The  change 
in  the  low-altitude  fields  is  almost  negligible  across  the  entire  band,  indicating  that 
the  results  are  accurate.  The  high-altitude  fields  differ  somewhat  above  20  kHz,  but 
the  change  is  still  small.  These  results  indicate  that  high-latitude  fields  are  correctly 
computed  with  1-km  grid  spacing. 
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Figure  3.  Test  of  FDTD  model  convergence  from  low-latitude  NPM  transmitter. 

The  panels  show  the  spectrum  of  low-altitude  (top)  and  high-altitude  (bottom) 
electric  fields  as  the  grid  spacing  in  the  simulation  is  reduced  by  25%.  The  change 
in  the  low-altitude  fields  is  almost  negligible  across  the  entire  band,  indicating  that 
the  results  are  accurate.  However,  the  high-altitude  fields  differ  somewhat  above  5 
kHz  and  differ  enormously  above  15  kHz. 

4.3.  High-Latitude  vs.  Low-Latitude  Model  Convergence — Complete 
Analysis 

To  quantify  the  connection  between  grid  spacing,  geomagnetic  latitude,  and 
simulation  accuracy,  we  ran  a  series  of  simulations  using  parameters  corresponding  to  the 
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NML  (high  latitude)  and  NPM  (low  latitude)  transmitters  with  grid  spacings  of  1000, 

500,  and  250  meters.  Finite  difference  simulations  are  guaranteed  to  eventually  converge 
to  the  correct  answer  as  this  grid  spacing  is  reduced,  but  it  is  not  always  easy  to  know 
what  grid  spacing  is  required  to  maintain  errors  below  a  specified  level.  We  have  found 
that,  for  a  given  uniform  grid  spacing,  the  low-altitude  fields  can  be  computed  correctly 
while  the  high-altitude  fields  are  dramatically  incorrect.  This  shows  that  a  nonuniform 
grid  approach,  in  which  low  altitudes  are  resolved  coarsely  and  only  the  highest  altitudes 
are  resolved  finely,  is  the  key  to  achieving  efficient  simulations  of  high-altitude  fields. 

At  high  latitudes  (i.e.,  NML),  500-m  grid  spacing  results  in  fairly  accurate  high-altitude 
fields  in  the  20  -  25-  Hz  band.  But  at  low  latitudes  (i.e.,  NPM),  250  m  and  smaller  grid 
spacing  is  required  to  achieve  the  same  accuracy.  The  convergence  in  all  cases  shows 
second-order  accuracy,  which  is  expected.  This  also  means  that  correct  answers  can  be 
estimated  by  extrapolating  the  convergence  to  its  limit. 

All  of  the  figures  below  show  the  amplitude  spectrum  of  one  electric  field 
component  (the  results  are  essentially  identical  for  all  E  and  H  components)  observed  at 
200-km  range  from  the  transmitter  (short  range  is  required  for  the  simulations  to  be 
manageable  for  fine  grids)  and  either  90-km  (low)  or  130-km  (high)  altitude.  Figure  4 
shows  the  low-altitude  (90  km)  field  spectrum  for  a  background  magnetic  field 
corresponding  to  the  NML  transmitter  (high  latitude)  for  3  different  grid  spacings.  The 
low-altitude  fields  are  almost  invariant  to  grid  spacing,  indicating  that  the  solution  has 
almost  converged,  even  at  1-km  spacing,  and,  consequently,  that  these  results  are  all 
numerically  correct.  If  we  treat  the  250-meter  solution  as  “correct,”  the  bottom  panel 
shows  the  relative  error  for  1-km  and  0.5-km  grids.  These  relative  errors  are  all  smaller 
than  10%  except  above  roughly  30  kHz. 
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Figure  4.  The  low-altitude  (90  km)  Held  spectrum  for  a  magnetic  latitude 
corresponding  to  the  NML  transmitter  (high  latitude)  for  3  different  grid  spacings. 


Figure  5  shows  the  comparison  for  die  same  field  component  at  130-km  altitude. 
At  low  frequencies  (<10  kHz),  the  solution  is  nearly  converged  for  1-km  grid  spacing. 

But  at  higher  frequencies,  particularly  the  target  frequencies  of  20-  to  30  kHz,  a  1-km 
grid  is  not  quite  enough  to  give  the  right  answer,  as  seen  by  the  almost  50%  change  in  the 
answer  when  the  grid  size  is  reduced  from  1  km  to  0.5  km.  The  step  from  0.5  to  0.25  km, 
however,  produces  only  a  small  change  (<10%)  in  the  solution,  indicating  that  by  0.5  km, 
the  solution  is  almost  correct  even  at  high  altitudes. 

Figure  5  (and  the  other  figures  below)  also  shows  clearly  that  the  solution  is 
converging  with  second-order  convergence,  in  which  the  error  shrinks  in  proportion  to 
the  square  of  the  grid  size.  In  other  words,  halving  the  grid  spacing  results  in  a  4x  error 
reduction  (see,  for  example,  how  the  error  drops  from  40%  to  10%  at  20  kHz  when  the 
grid  goes  from  1  km  to  0.5  km).  This  is  as  it  should  be,  but  it  is  always  nice  to  get 
confirmation. 

Importantly,  because  this  convergence  is  well-understood  mathematically,  it  can 
be  exploited  to  estimate  the  “correct”  answer  from  simulations  that  have  not  quite 
converged  to  the  correct  answer.  This  idea  is  called  the  “deferred  approach  to  the  limit” 
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and  is  used  frequently  in  a  variety  of  numerical  calculations.  For  example,  in  the  upper 
right  panel,  we  can  see  that  the  fields  are  increasing  steadily  as  the  grid  size  is  reduced, 
and  we  know  that  the  error  is  decreasing  by  a  factor  of  4  for  every  halving  of  the  grid 
size.  Knowing  this,  the  limit  of  this  convergence  can  be  estimated  from  simulations  at 
only  two  grid  sizes,  each  of  which  has  not  yet  converged.  In  this  way,  the  desired  result 
can  be  obtained  using  relatively  coarse  and,  therefore,  much  faster  computations. 
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Figure  5.  The  high-altitude  (130  km)  field  spectrum  for  a  magnetic  latitude 
corresponding  to  the  NML  transmitter  (high  latitude)  for  3  different  grid  spacings. 

The  story  is  somewhat  different  for  low-latitude  transmitters,  as  we  have  already 
found  and  reported.  Figure  6  shows,  in  the  same  format  as  Figure  4,  the  low-altitude  (90 
km)  fields  produced  by  the  NPM  transmitter.  For  the  same  sequence  of  grid  spacings,  the 
solution  is  not  nearly  as  converged  as  was  found  for  the  high-latitude  transmitter.  A  500- 
meter  grid  is  needed  to  reach  the  same  relative  error  (1-2%)  that  a  1-km  grid  yielded  for 
NML.  Nevertheless,  reasonably  correct  answers  are  obtained  even  for  the  larger  grid 
spacings  considered  here.  Note,  again,  the  second-order  convergence — for  example,  the 
error  drops  from  8%  to  2%  at  about  17  kHz  when  the  grid  size  is  halved. 
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Figure  6.  The  low-altitude  (90  km)  field  spectrum  for  a  magnetic  latitude 
corresponding  to  the  NPM  transmitter  (low  latitude)  for  3  different  grid  spacings. 


The  convergence  of  the  high-altitude  fields  is  especially  striking,  as  shown  in 
Figure  7.  Even  though  a  1-km  grid  yields  the  low-altitude  (90  km)  fields  with  reasonable 
accuracy  up  to  30  kHz,  the  high-altitude  (130  km)  fields  are  grossly  underestimated 
above  about  5  kHz  for  the  same  grid.  The  solution  even  changes  significantly  (about  a 
factor  of  2)  when  the  grid  size  is  reduced  from  500  to  250  meters,  indicating  that  an  even 
finer  grid  is  required  for  an  answer  close  to  convergence.  We  expect  that  the  deferred 
approach  to  the  limit  described  above  may  be  especially  valuable  for  these  low-latitude 
cases  in  which  the  required  grid  spacing  is  fine. 

Together,  these  results  quantify  rather  precisely  the  relationship  between  grid  size 
and  error  at  both  low-  and  high-latitudes  and  low-  and  high-altitudes.  Importantly,  they 
confirm  that  the  computations  that  we  have  already  reported  for  the  high-latitude  NML 
transmitter,  which  used  a  1-km  grid,  are  reasonably  accurate  but  perhaps  underestimate 
the  high-altitude  fields  by  roughly  30%. 
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Figure  7.  The  high-altitude  (130  km)  field  spectrum  for  a  magnetic  latitude 
corresponding  to  the  NPM  transmitter  (low  latitude)  for  3  different  grid  spacings 


4.4.  Time  Harmonic  Code  and  Initial  Results 

As  mentioned  above,  we  are  developing  a  time  harmonic  finite  difference  code  to 
complement  the  time  domain  code  that  has  produced  the  already-reported  results  from 
this  project.  A  time  harmonic  code  will  have  several  advantages.  For  single  frequency 
computations,  it  should  be  significantly  faster  than  the  time  domain  code.  It  also  enables 
variable  grid  resolution  that  can  be  implemented  easily,  which  should  also  increase  model 
efficiency  by  letting  fine  grid  spacing  be  used  only  where  it  is  needed  (high  altitude,  in 
general). 

Here,  we  report  results  comparing  frequency  domain  and  time  domain  simulations 
for  high-latitude  locations.  The  plots  that  follow  compare  low-altitude  (0  km)  and  high- 
altitude  (130  km)  VLF  electric  fields  computed  by  the  time  domain  code  and  the  new 
frequency  domain  code.  The  frequency  domain  code  is  functional  although  it  is  still  very 
much  in  development.  Because  of  limitations  of  the  Matlab  matrix  solvers,  we  use,  in 
this  preliminary  version,  the  horizontal  range  of  the  simulations  which  is  limited  to  about 


12 


100  km.  Obviously,  this  is  much  shorter  than  we  ultimately  need,  but  we  emphasize  that 
this  is  purely  a  limitation  on  the  solvers,  not  on  the  technique  or  even  the  code. 
Transitioning  this  code  to  other  computers  and  solvers  should  be  a  very  straightforward 
process.  The  simulations  shown  below  are  also  for  high-latitude  sources  for  simplicity. 

Figure  8  shows  the  vertical  electric  field  vs.  range  at  ground  level  for  a  10-kHz 
source  (radiated  power  is  not  important  as  long  as  the  sources  are  identical  in  both 
simulations)  .  The  plot  shows  field  amplitude  for  both  frequency  domain  (FD)  and  time 
domain  codes  (TD)  computed  with  two  different  levels  of  spatial  discretization  (500  m 
and  1000  m).  Aside  from  some  small  discrepancies  within  10  km  of  the  source,  the  two 
are  in  excellent  agreement.  This  also  confirms  that  low-altitude  fields  are  essentially 
converged  at  1-km  discretization. 
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Figure  8.  Low-altitude  fields  produced  by  a  10-kHz  source  computed  using  the  FD 
and  TD  codes.  The  agreement  is  excellent,  validating  the  new  FD  code. 

Figure  9  shows  all  three  electric  field  components  at  high  altitude  (130  km)  for 
the  same  source  computed  with  both  codes.  The  FD  code  results  change  somewhat  when 
the  spatial  discretization  is  reduced  to  500  m,  indicating  that  the  simulation  is  not  quite 
converged  at  1000-m  discretization.  However,  the  500-m  TD  and  FD  runs  are  again  in 
excellent  quantitative  agreement,  confirming  the  correctness  of  the  FD  simulation  results. 
The  component  of  E  parallel  to  the  background  magnetic  field  (Ez,  in  this  case)  is  much 
smaller  than  the  other  components,  as  expected.  The  two  horizontal  components  have 
equal  amplitudes  (shown)  and  are  in  phase  quadrature,  showing  that  the  upward  traveling 
fields  are  circularly  polarized,  also  as  expected. 

Figures  10  and  1 1  show  the  same  comparisons  except  for  a  20-kHz  source 
frequency.  The  low-altitude  fields  agree  well,  which  is  consistent  with  the  TD  code 
convergence  analysis.  The  high-altitude  fields  are  generally  in  good  quantitative 
agreement,  but  not  to  the  same  degree  as  the  10-kHz  fields.  The  FD  results  also  show  a 
major  change  in  the  shape  of  the  field  profile  from  1000-m  to  500-m  discretization.  A 
further  decrease  in  the  FD  discretization  might  be  required  for  more  quantitatively 
accurate  results.  The  most  important  quantity,  the  general  field  amplitude,  agrees  well 
between  the  TD  and  FD  code. 
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Figure  9.  High-altitude  fields  produced  by  a  10-kHz  source  computed  using  the  FD 
and  TD  codes.  The  agreement  is  again  excellent. 


E  amplitude  at  ground  (f = 20  kHz) 


Figure  10.  Low-altitude  fields  produced  by  a  20-kHz  source  computed  using  the  FD 
and  TD  codes. 
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Figure  11.  High-altitude  fields  produced  by  a  20-kHz  source  computed  using  the 
FD  and  TD  codes. 


The  important  difference  between  these  codes  is  their  speed,  which  is  shown  in 
Table  1.  The  TD  code  computes  all  frequencies  simultaneously,  which  has  some 
advantages  for  broadband  fields  but  may  not  be  maximally  efficient  if  only  one  frequency 
or  several  frequencies  are  desired.  The  FD  code  computes  the  field  for  only  one  specific 
frequency,  and  is  consequently  faster.  The  two  codes  were  run  on  comparable  machines. 
Each  code’s  absolute  run  time  could  be  improved  by  running  on  faster  platforms,  but  the 
ratio  between  the  two  should  be  relatively  invariant. 


Table  1.  Comparison  of  TD  and  FD  Code  Run  Times.  The  FD  Code  is 
Approximately  90  Times  Faster 


Code  Type 

Spatial  Discretization 

Run  Time 

TD 

1000  m 

1332  s 

FD 

1000  m 

16  s 

TD 

500  m 

9684  s 

FD 

500  m 

129  s 
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As  Table  1  shows,  the  ED  code  is  approximately  90  times  faster  than  the  TD  code 
for  both  discretizations.  This  is  a  major  improvement  and  justifies  our  further  effort  in 
developing  the  FD  code.  We  also  emphasize  that  the  FD  code  makes  multiple  grid  scales 
relatively  straightforward.  This  will  also  improve  code  efficiency  by  enabling  us  to  place 
spatial  resolution  only  where  it  is  needed. 

5.  CONCLUSIONS 

Despite  challenges  related  to  the  arrival  of  funds  to  support  this  project, 
significant  technical  progress  has  been  made.  A  complete  end-to-end  run  of  the  high- 
altitude  VLF  power  prediction  model  has  been  completed  using  our  modeled  130-km 
altitude  VLF  power  predictions  as  the  input  to  the  higher  altitude  ray-tracing  code.  This 
resulted  in  good  agreement  with  high-altitude  VLF  field  measurements  from  the  IMAGE 
satellite.  Previous  end-to-end  runs  like  this  which  used  semi-empirical  ionospheric 
penetration  had  revealed  a  ~20-dB  discrepancy  between  the  measured  and  predicted 
fields.  Our  model  does  not  require  any  such  correction  factor  and  appears  to  show  that 
this  discrepancy  originates  in  the  ionospheric  penetration  calculation. 

We  also  better  understand  the  convergence  and  accuracy  issues  associated  with 
the  magnetic  latitude  of  the  source.  A  series  of  simulations  using  parameters 
corresponding  to  the  NML  (high  latitude)  and  NPM  (low  latitude)  transmitters  with  grid 
spacings  of  1000,  500,  and  250  meters  showed  that,  for  a  given  uniform  grid  spacing,  the 
low-altitude  fields  can  be  computed  correctly  while  the  high-altitude  fields  are 
dramatically  incorrect.  This  shows  that  a  nonuniform  grid  approach,  in  which  low 
altitudes  are  resolved  coarsely  and  only  the  highest  altitudes  are  resolved  finely,  is  the 
key  to  achieving  efficient  simulations  of  high-altitude  fields.  At  high  latitudes  (i.e., 
NML),  500-m  grid  spacing  results  in  fairly  accurate  high-altitude  fields  in  the  20-25- 
kHz  band.  But,  at  low  latitudes  (i.e.,  NPM),  250  m  and  smaller  grid  spacing  is  required 
to  achieve  the  same  accuracy.  The  convergence  in  all  cases  shows  second-order  accuracy, 
which  is  expected.  This  also  means  that  correct  answers  can  be  estimated  by 
extrapolating  the  convergence  to  its  limit. 

Lastly,  the  realization  that  very  fine  grid  spacings  are  required  to  achieve  accurate 
results  at  low  latitudes  has  led  us  to  develop  a  more  efficient  frequency  domain  (i.e., 
single  frequency)  VLF  solver.  Comparisons  with  results  from  the  time  domain  code 
show  that  our  preliminary  frequency  domain  code  gives  answers  in  quantitative 
agreement  with  the  time  domain  code  but  is  about  90  times  faster  for  single-frequency 
computations.  This  major  enhancement  in  computational  speed  justifies  the  further 
development  needed  for  this  solver. 

The  ongoing  code  development  will  enable  accurate  predictions  of  the 
ionospheric  wave  fields  and  their  statistical  uncertainty.  The  end  result  will  be  the 
capability  to  predict  more  accurately  the  injection  of  VLF  wave  energy  into  the 
magnetosphere,  which  is  a  key  component  in  the  accurate  modeling  of  the  influence  of 
VLF  on  natural  radiation  belt  dynamics. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AFRL 

FDTD 

LWPC 

NPML 

VLF 


Air  Force  Research  Laboratory 
Finite  Difference  Time  Domain 
Long  Wave  Propagation  Capability 
Nearly  Perfect  Matched  Layer 
Very  Low  Frequency 
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